Load data

vcs <- rio::import("data_complete_all_anonym_allZ_unfiltered.xlsx")
var_label(vcs$dominance) <- "Dominance"
var_label(vcs$neuro) <- "Neuroticism"
var_label(vcs$agree) <- "Agreeableness"
var_label(vcs$extra) <- "Extraversion"
var_label(vcs$openn) <- "Openness"
var_label(vcs$consc) <- "Conscientiousness"
var_label(vcs$soi_full) <- "Unrestricted sociosexuality"
var_label(vcs$f0) <- "Voice pitch"
var_label(vcs$pf) <- "Formants"
vcs$sex_c <- vcs$sex
vcs$sex <- factor(if_else(vcs$sex == 1, "male", "female"))
contrasts(vcs$sex) <- contr.helmert(2)
var_label(vcs$sex) <- "Sex"
set.seed(1)
var_label(vcs$age) <- "Age"
vcs <- vcs %>% 
  mutate(
    age = if_else(dataset == 9, 20, age)/10,
    age_se = if_else(dataset == 9, 3, 0.5)/10
  )
vcs <- vcs %>% filter(!is.na(f0), !is.na(pf), !is.na(age))

warmup <- 2000
iter <- warmup + 2000
chains <- 4
control <- list(adapt_delta = 0.99)
priors <- c(
  prior(normal(0, 3), class = b)
)

variable_labels <- c("Intercept"= "Intercept",  "sex1" = "Sex [male]", "bsp_meageage_se" = "Age±SE", "f0" = "Voice pitch (f0)", "pf" = " Formant position (Pf)", "age" = "Age", "f0:sex1" = "Voice pitch (f0) × Sex")
effect_labels <- c("b_Intercept"= "Intercept",  "b_sex1" = "Sex [male]", "bsp_meageage_se" = "Age±SE", "b_f0" = "Voice pitch (f0)", "b_pf" = " Formant position (Pf)", "b_age" = "Age", "b_f0.sex1" = "Voice pitch (f0) × Sex")

Hypothesis 1: Dominance

Participants with lower voice pitch will have a more dominant personality.

Interpretation: The data are are consistent with substantial linear negative relationships between f0 and dominance. The 89% HDI for f0 falls entirely outside the ROPE, but the HDI for pf falls almost entirely within it. This evidence is consistent with a non-negligible association, where people with deeper voices are more dominant (after adjusting for age and gender). Further research is needed to verify if the association with pf is truly negligible in size.

Visually diagnose non-linearity

ggplot(vcs, aes(f0, dominance)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(pf, dominance)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(age*10, dominance)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Decisions upon visual diagnosis: Visual diagnostic show no clear sign of nonlinearities. We will fit linear effects for f0, pf, and age. It is not necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is not included in these analyses. Visually, there could be an interaction between sex and f0.

Models

h1_simple <- brm(dominance ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/dominance/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 


h1_after_diagnosis <- brm(dominance ~ f0*sex + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, cores = chains, chains = chains,
          prior = priors, 
          control = control, save_mevars = TRUE,
          file = "models/dominance/h1_nonlinear") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
# conditional_smooths(h1_after_diagnosis)

compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
##                    elpd_diff se_diff
## h1_simple           0.0       0.0   
## h1_after_diagnosis -1.0       0.3
h1_sex_mod_dataset_varying <- brm(dominance ~ f0 + sex + age + 
                                    (1 + f0 + sex || dataset), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/dominance/h1_sex_mod_dataset_varying2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_after_diagnosis_melsm <- brm(
  bf(dominance ~ f0 + sex + age + (1 | dataset),
     sigma ~ f0 + sex + age + (1 | dataset)), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/dominance/h1_after_diagnosis_melsm") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
  

compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying,
                      h1_after_diagnosis_melsm)
compare_models
##                            elpd_diff se_diff
## h1_after_diagnosis_melsm    0.0       0.0   
## h1_simple                  -3.9       3.9   
## h1_sex_mod_dataset_varying -4.4       4.0   
## h1_after_diagnosis         -4.9       4.0
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying,
                      h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>% 
  map(~ { x <- as_tibble(bayes_R2(.))
          x$loo_R2 <- loo_R2(.)
          x
          }) %>%
  bind_rows(.id="model") %>% 
  arrange(desc(loo_R2))

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  Dominance
Predictors Estimates CI (95%)
Intercept -0.07 -0.75;0.60
Voice pitch (f0) -0.27 -0.45;-0.09
Formant position (Pf) -0.02 -0.20;0.16
Sex [male] -0.31 -0.53;-0.09
Age 0.04 -0.08;0.16
Random Effects
σ2 0.05
Ï„00 0.96
ICC 0.05
N dataset 4
Observations 985
Marginal R2 / Conditional R2 0.069 / 0.069 
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter H0 inside ROPE 89% HDI
Intercept Undecided 31% [-0.59; 0.40]
Voice pitch (f0) Rejected 0% [-0.41;-0.11]
Formant position (Pf) Undecided 79% [-0.16; 0.13]
Sex [male] Rejected 0% [-0.48;-0.12]
Age Undecided 89% [-0.06; 0.13] 
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Model details

summary(h1_simple)
## Warning: There were 3 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: dominance ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 985) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.51      0.36     0.16     1.52 1.00     2354     2909
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.07      0.33    -0.75     0.60 1.00     2938     3185
## f0           -0.27      0.09    -0.45    -0.09 1.00     5462     5238
## pf           -0.02      0.09    -0.20     0.16 1.00     6544     5231
## sex1         -0.31      0.11    -0.53    -0.09 1.00     4937     4930
## age           0.04      0.06    -0.08     0.16 1.00     6745     5282
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.93     1.01 1.00     7476     5083
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis)
## Warning: There were 7 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: dominance ~ f0 * sex + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 985) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.51      0.41     0.16     1.58 1.00     2139     2550
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.05      0.36    -0.79     0.63 1.00     3316     3248
## f0           -0.26      0.10    -0.45    -0.07 1.00     6176     5797
## sex1         -0.30      0.12    -0.52    -0.07 1.00     5648     5260
## pf           -0.02      0.09    -0.20     0.16 1.00     7408     5498
## age           0.04      0.06    -0.09     0.16 1.00     7674     5247
## f0:sex1       0.02      0.10    -0.17     0.21 1.00     8212     5641
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.93     1.02 1.00     8345     5503
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(h1_after_diagnosis_melsm)

# conditional_smooths(h1_after_diagnosis)
summary(h1_sex_mod_dataset_varying)
## Warning: There were 21 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: dominance ~ f0 + sex + age + (1 + f0 + sex || dataset) 
##    Data: vcs (Number of observations: 985) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.60      0.46     0.18     1.82 1.00     4020     5829
## sd(f0)            0.18      0.24     0.00     0.80 1.00     4385     4587
## sd(sex1)          0.31      0.41     0.01     1.48 1.00     2541     2752
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.09      0.41    -0.94     0.74 1.00     4552     2993
## f0           -0.26      0.17    -0.57     0.05 1.00     6829     4843
## sex1         -0.28      0.27    -0.78     0.29 1.00     3858     2606
## age           0.04      0.06    -0.08     0.16 1.00    12997     8476
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.93     1.01 1.00    12338     5464
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis_melsm)
## Warning: There were 4 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: dominance ~ f0 + sex + age + (1 | dataset) 
##          sigma ~ f0 + sex + age + (1 | dataset)
##    Data: vcs (Number of observations: 985) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.52      0.39     0.17     1.59 1.00     3578     4584
## sd(sigma_Intercept)     0.14      0.17     0.01     0.56 1.00     2581     3435
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          -0.12      0.35    -0.78     0.53 1.00     5000     4920
## sigma_Intercept     0.05      0.16    -0.28     0.34 1.00     4570     4376
## f0                 -0.30      0.10    -0.49    -0.11 1.00     8536     8368
## sex1               -0.33      0.10    -0.51    -0.14 1.00     8776     7747
## age                 0.05      0.06    -0.07     0.17 1.00    11548     8335
## sigma_f0            0.11      0.07    -0.03     0.24 1.00     8723     7404
## sigma_sex1          0.11      0.07    -0.03     0.25 1.00     8507     7525
## sigma_age          -0.02      0.05    -0.12     0.07 1.00     7026     7494
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
Explained variance, regular and leave-one-out.
model Estimate Est.Error Q2.5 Q97.5 loo_R2
h1_after_diagnosis_melsm 0.0704834 0.0144778 0.0435837 0.1001937 0.0545025
h1_simple 0.0699320 0.0150681 0.0425830 0.1009986 0.0530509
h1_sex_mod_dataset_varying 0.0715111 0.0148753 0.0439116 0.1018806 0.0519595
h1_after_diagnosis 0.0707440 0.0151169 0.0430927 0.1016319 0.0510680 
kable(coef(h1_sex_mod_dataset_varying)$dataset[, ,"f0"], caption = "Estimated association between f0 and outcome by dataset")
Estimated association between f0 and outcome by dataset
Estimate Est.Error Q2.5 Q97.5
2 -0.2973410 0.1170703 -0.5362739 -0.0785982
3 -0.2311627 0.1267894 -0.4659943 0.0393573
7 -0.2364938 0.1564241 -0.5189253 0.1130254
8 -0.2943222 0.1467621 -0.6110407 -0.0168347

Hypothesis 2: Extraversion

Participants with lower voice pitch will score higher on extraversion

Interpretation: The data are are consistent with a substantial linear negative relationship between f0 and extraversion. The 89% HDI for f0 excludes and falls entirely outside the ROPE. The 89% of the HDI for pf overlaps zero and the ROPE boundaries. This evidence is consistent with an association, where people with deeper voices have higher extraversion (after adjusting for age and gender), but further research is needed to test whether the association with pf is negligible in size.

Visually diagnose non-linearity

ggplot(vcs, aes(f0, extra)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour =  "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(pf, extra)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(age*10, extra)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Decisions upon visual diagnosis: Visual diagnostic shows signs of a potential interaction by sex coupled with nonlinearity for pf. We will fit separate splines by sex for pf. It is necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is included in these analyses.

Models

h1_simple <- brm(extra ~ f0 + pf + sex + me(age, age_se) + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/extra/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 


h1_after_diagnosis <- brm(extra ~ f0 + s(pf, by = sex) + sex + me(age, age_se) + (1 | dataset), data = vcs,
          iter = iter + 1000, warmup = warmup + 1000, cores = chains, chains = chains,
          prior = priors, 
          control = control, save_mevars = TRUE,
          file = "models/extra/h1_nonlinear2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
## Warning: Found 1 observations with a pareto_k > 0.7 in model '.'. It is recommended to set 'reloo = TRUE' in
## order to calculate the ELPD without the assumption that these observations are negligible. This will refit the
## model 1 times to compute the ELPDs for the problematic observations directly.
conditional_smooths(h1_after_diagnosis)

compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
##                    elpd_diff se_diff
## h1_simple           0.0       0.0   
## h1_after_diagnosis -1.4       1.3
h1_sex_mod_dataset_varying <- brm(extra ~ f0 + sex + me(age, age_se) + 
                                    (1 + f0 + sex || dataset), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/extra/h1_sex_mod_dataset_varying2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_after_diagnosis_melsm <- brm(
  bf(extra ~ f0 + sex + me(age, age_se) + (1 | dataset),
     sigma ~ f0 + sex + me(age, age_se) + (1 | dataset)), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/extra/h1_after_diagnosis_melsm") %>%
  add_criterion("loo") %>%
  add_criterion("bayes_R2") %>%
  add_criterion("loo_R2")


compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying ,h1_after_diagnosis_melsm
                      )
compare_models
##                            elpd_diff se_diff
## h1_after_diagnosis_melsm     0.0       0.0  
## h1_simple                  -14.7       5.3  
## h1_sex_mod_dataset_varying -15.1       5.4  
## h1_after_diagnosis         -16.1       5.5
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying,h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>% 
  map(~ { x <- as_tibble(bayes_R2(.))
          x$loo_R2 <- loo_R2(.)
          x
          }) %>%
  bind_rows(.id="model") %>% 
  arrange(desc(loo_R2))

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  Extraversion
Predictors Estimates CI (95%)
Intercept 0.17 -0.22;0.57
Voice pitch (f0) -0.23 -0.38;-0.08
Formant position (Pf) 0.08 -0.05;0.21
Sex [male] -0.30 -0.48;-0.11
Age -0.10 -0.22;0.01
Random Effects
σ2 0.03
Ï„00 0.97
ICC 0.03
N dataset 7
Observations 1433
Marginal R2 / Conditional R2 0.055 / 0.055 
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter H0 inside ROPE 89% HDI
Intercept Undecided 31% [-0.14; 0.48]
Voice pitch (f0) Rejected 0% [-0.35;-0.11]
Formant position (Pf) Undecided 64% [-0.03; 0.19]
Sex [male] Rejected 0% [-0.45;-0.15]
Age Undecided 46% [-0.20;-0.02] 
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Model details

summary(h1_simple)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: extra ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1433) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.35      0.15     0.16     0.74 1.00     2181     3466
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.17      0.20    -0.22     0.57 1.00     3364     4103
## f0           -0.23      0.08    -0.38    -0.08 1.00     5308     5343
## pf            0.08      0.07    -0.05     0.21 1.00     5653     5187
## sex1         -0.30      0.09    -0.48    -0.11 1.00     4663     5056
## age          -0.10      0.06    -0.22     0.01 1.00     6253     5299
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.94     1.01 1.00     7346     5546
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: extra ~ f0 + s(pf, by = sex) + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 1433) 
## Samples: 4 chains, each with iter = 5000; warmup = 3000; thin = 1;
##          total post-warmup samples = 8000
## 
## Smooth Terms: 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(spfsexfemale_1)     0.78      0.71     0.03     2.63 1.00     5276     4177
## sds(spfsexmale_1)       1.00      0.85     0.04     3.10 1.00     4908     4298
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.34      0.15     0.16     0.74 1.00     2973     4582
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           0.15      0.20    -0.25     0.55 1.00     3408     3909
## f0                 -0.23      0.08    -0.38    -0.08 1.00     7969     6731
## sex1               -0.31      0.10    -0.50    -0.13 1.00     6763     5319
## spf:sexfemale_1     1.01      1.70    -2.53     4.42 1.00     6309     5392
## spf:sexmale_1       0.68      1.73    -2.39     4.60 1.00     6668     5360
## meageage_se        -0.11      0.06    -0.22     0.00 1.00     8086     6273
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.94     1.01 1.00    13235     5381
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(h1_after_diagnosis_melsm)

# conditional_smooths(h1_after_diagnosis_melsm)
summary(h1_sex_mod_dataset_varying)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: extra ~ f0 + sex + me(age, age_se) + (1 + f0 + sex || dataset) 
##    Data: vcs (Number of observations: 1433) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.37      0.17     0.17     0.80 1.00     4339     6085
## sd(f0)            0.09      0.09     0.00     0.31 1.00     5237     6722
## sd(sex1)          0.11      0.12     0.00     0.42 1.00     4314     5680
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.18      0.22    -0.24     0.61 1.00     5476     5979
## f0             -0.23      0.09    -0.41    -0.05 1.00     9479     7800
## sex1           -0.32      0.11    -0.52    -0.09 1.00     8368     7560
## meageage_se    -0.10      0.06    -0.22     0.01 1.00    12329     9326
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.02     0.94     1.01 1.00    21123     8908
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis_melsm)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: extra ~ f0 + sex + me(age, age_se) + (1 | dataset) 
##          sigma ~ f0 + sex + me(age, age_se) + (1 | dataset)
##    Data: vcs (Number of observations: 1433) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.36      0.15     0.17     0.76 1.00     4273     6355
## sd(sigma_Intercept)     0.14      0.07     0.05     0.31 1.00     3801     6797
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             0.17      0.20    -0.21     0.57 1.00     4649     5661
## sigma_Intercept       0.04      0.12    -0.21     0.28 1.00     7281     8031
## f0                   -0.20      0.08    -0.35    -0.05 1.00    10353     9243
## sex1                 -0.32      0.08    -0.47    -0.17 1.00    10062     8618
## sigma_f0              0.09      0.06    -0.02     0.20 1.00    10765     8576
## sigma_sex1            0.12      0.06     0.00     0.23 1.00     9947     8834
## meageage_se          -0.10      0.05    -0.20     0.00 1.00    11772     8801
## sigma_meageage_se    -0.01      0.04    -0.10     0.08 1.00     8560     8003
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
Explained variance, regular and leave-one-out.
model Estimate Est.Error Q2.5 Q97.5 loo_R2
h1_after_diagnosis_melsm 0.0525347 0.0113320 0.0318093 0.0762478 0.0415888
h1_simple 0.0550861 0.0112503 0.0345326 0.0781500 0.0408287
h1_sex_mod_dataset_varying 0.0569871 0.0114281 0.0360700 0.0808472 0.0400963
h1_after_diagnosis 0.0575437 0.0113746 0.0364921 0.0815037 0.0386056 
kable(coef(h1_sex_mod_dataset_varying)$dataset[, ,"f0"], caption = "Estimated association between f0 and outcome by dataset")
Estimated association between f0 and outcome by dataset
Estimate Est.Error Q2.5 Q97.5
2 -0.1996990 0.0918685 -0.3709783 -0.0125423
3 -0.1966691 0.0985113 -0.3871866 0.0008238
4 -0.2098008 0.1069189 -0.4139233 0.0115944
7 -0.2139830 0.1162529 -0.4393242 0.0251057
8 -0.2679709 0.1263647 -0.5629834 -0.0522468
9 -0.2414606 0.1034007 -0.4538604 -0.0427239
10 -0.2546173 0.1271442 -0.5420006 -0.0283154 
equivalence_test(h1_after_diagnosis_melsm, range = c(-0.1, 0.1))
## # Test for Practical Equivalence
## 
##   ROPE: [-0.10 0.10]
## 
## Parameter         |        H0 | inside ROPE |       89% HDI
## -----------------------------------------------------------
## Intercept         | Undecided |     31.44 % | [-0.14  0.47]
## sigma_Intercept   | Undecided |     62.09 % | [-0.14  0.25]
## f0                | Undecided |      3.99 % | [-0.32 -0.08]
## sex1              |  Rejected |      0.00 % | [-0.45 -0.20]
## sigma_f0          | Undecided |     53.07 % | [ 0.01  0.19]
## sigma_sex1        | Undecided |     38.34 % | [ 0.02  0.21]
## meageage_se       | Undecided |     50.59 % | [-0.18 -0.02]
## sigma_meageage_se |  Accepted |    100.00 % | [-0.08  0.06]

Hypothesis 3: Agreeableness

Participants with lower voice pitch will score higher on agreeableness.

Interpretation: The credible intervals for f0 and pf include zero, but the 89% HDI does not fall entirely within the ROPE. This evidence is consistent with a negligible association between f0/pf and agreeableness. Further research is needed to verify if the associations are truly negligible in size.

Visually diagnose non-linearity

ggplot(vcs, aes(f0, agree)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour =  "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(pf, agree)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(age*10, agree)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Decisions upon visual diagnosis: Visual diagnostic shows signs of nonlinearity for f0 and pf, but no clear interactions by sex. We will fit splines for f0, and pf. It is necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is included in these analyses.

Models

h1_simple <- brm(agree ~ f0 + pf + sex + me(age, age_se) + (1 | dataset), data = vcs,
          iter = iter + 1000, warmup = warmup + 1000, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/agree/h1_simple2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 



h1_after_diagnosis <- brm(agree ~ f0 * sex + pf + sex + me(age, age_se) + (1 | dataset), data = vcs,
                          iter = iter + 1500, warmup = warmup + 1500, cores = chains, chains = chains,
                          prior = priors, 
                          control = control, save_mevars = TRUE,
                          file = "models/agree/h1_nonlinear4") %>% 
    add_criterion("loo") %>% 
    add_criterion("bayes_R2") %>% 
    add_criterion("loo_R2") 

compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
##                    elpd_diff se_diff
## h1_after_diagnosis  0.0       0.0   
## h1_simple          -2.7       2.6
h1_sex_mod_dataset_varying <- brm(agree ~ f0 * sex + sex + me(age, age_se) + 
                                    (1 + f0 + sex || dataset), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/agree/h1_sex_mod_dataset_varying2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_after_diagnosis_melsm <- brm(
  bf(agree ~ f0 * sex + sex + me(age, age_se) + (1 | dataset),
     sigma ~ f0 * sex + sex + me(age, age_se) + (1 | dataset)), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/agree/h1_after_diagnosis_melsm2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
  

compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying,
                      h1_after_diagnosis_melsm)
compare_models
##                            elpd_diff se_diff
## h1_after_diagnosis_melsm     0.0       0.0  
## h1_sex_mod_dataset_varying  -4.4       5.8  
## h1_after_diagnosis          -7.4       4.8  
## h1_simple                  -10.1       5.4
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying,
                      h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>% 
  map(~ { x <- as_tibble(bayes_R2(.))
          x$loo_R2 <- loo_R2(.)
          x
          }) %>%
  bind_rows(.id="model") %>% 
  arrange(desc(loo_R2))

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  Agreeableness
Predictors Estimates CI (95%)
Intercept 0.08 -0.38;0.56
Voice pitch (f0) 0.02 -0.13;0.17
Formant position (Pf) 0.04 -0.09;0.17
Sex [male] -0.01 -0.18;0.17
Age±SE -0.06 -0.18;0.05
Random Effects
σ2 0.11
Ï„00 0.90
ICC 0.11
N dataset 7
Observations 1434
Marginal R2 / Conditional R2 0.117 / 0.117 
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter H0 inside ROPE 89% HDI
Intercept Undecided 37% [-0.32; 0.43]
Voice pitch (f0) Undecided 91% [-0.10; 0.14]
Formant position (Pf) Undecided 85% [-0.06; 0.15]
Sex [male] Undecided 83% [-0.14; 0.15]
Age±SE Undecided 77% [-0.15; 0.03] 
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Model details

summary(h1_simple)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: agree ~ f0 + pf + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 5000; warmup = 3000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.47      0.20     0.24     0.98 1.00     2442     3519
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.08      0.24    -0.38     0.56 1.00     2346     3356
## f0              0.02      0.07    -0.13     0.17 1.00     7547     6037
## pf              0.04      0.06    -0.09     0.17 1.00     9301     6292
## sex1           -0.01      0.09    -0.18     0.17 1.00     6072     5617
## meageage_se    -0.06      0.06    -0.18     0.05 1.00     8287     6177
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.94      0.02     0.91     0.98 1.00    14245     6162
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: agree ~ f0 * sex + pf + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 5500; warmup = 3500; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.47      0.19     0.25     0.97 1.00     2345     3921
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.30      0.25    -0.19     0.80 1.00     3047     3641
## f0              0.08      0.08    -0.07     0.23 1.00     7448     6280
## sex1            0.09      0.10    -0.10     0.28 1.00     6250     5985
## pf              0.04      0.07    -0.08     0.17 1.00     9015     6260
## f0:sex1         0.20      0.07     0.06     0.35 1.00    11556     6633
## meageage_se    -0.08      0.06    -0.19     0.04 1.00     9122     5402
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.94      0.02     0.91     0.98 1.00    13050     5738
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(h1_after_diagnosis_melsm)

# conditional_smooths(h1_after_diagnosis_melsm)
summary(h1_sex_mod_dataset_varying)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: agree ~ f0 * sex + sex + me(age, age_se) + (1 + f0 + sex || dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.52      0.22     0.25     1.08 1.00     3931     5968
## sd(f0)            0.16      0.12     0.01     0.46 1.00     3585     4278
## sd(sex1)          0.13      0.14     0.00     0.50 1.00     3600     4159
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.24      0.27    -0.29     0.79 1.00     4163     5688
## f0              0.06      0.11    -0.16     0.27 1.00     6897     5709
## sex1            0.04      0.13    -0.22     0.28 1.00     5010     4462
## f0:sex1         0.14      0.08    -0.03     0.31 1.00     9624     8793
## meageage_se    -0.06      0.06    -0.17     0.05 1.00    11263     8694
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.94      0.02     0.90     0.97 1.00    16769     9043
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis_melsm)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: agree ~ f0 * sex + sex + me(age, age_se) + (1 | dataset) 
##          sigma ~ f0 * sex + sex + me(age, age_se) + (1 | dataset)
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.47      0.20     0.24     0.97 1.00     3897     6032
## sd(sigma_Intercept)     0.11      0.06     0.02     0.26 1.00     3596     3226
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             0.28      0.24    -0.18     0.76 1.00     4997     6522
## sigma_Intercept       0.16      0.13    -0.10     0.41 1.00     9090     8464
## f0                    0.09      0.08    -0.06     0.24 1.00     9620     8545
## sex1                  0.06      0.08    -0.11     0.23 1.00     8988     8427
## f0:sex1               0.18      0.07     0.03     0.32 1.00    13830     9184
## sigma_f0             -0.08      0.06    -0.20     0.05 1.00     8689     8631
## sigma_sex1           -0.08      0.07    -0.22     0.05 1.00     8625     8025
## sigma_f0:sex1        -0.02      0.06    -0.14     0.09 1.00    14660     9598
## meageage_se          -0.07      0.05    -0.17     0.02 1.00    15033     9269
## sigma_meageage_se    -0.09      0.04    -0.17    -0.01 1.00     9195     8811
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
Explained variance, regular and leave-one-out.
model Estimate Est.Error Q2.5 Q97.5 loo_R2
h1_sex_mod_dataset_varying 0.1280173 0.0154562 0.0982128 0.1586401 0.1109300
h1_after_diagnosis_melsm 0.1198456 0.0157054 0.0899097 0.1514868 0.1085033
h1_after_diagnosis 0.1220400 0.0151115 0.0931632 0.1525303 0.1074057
h1_simple 0.1174470 0.0150143 0.0892506 0.1477528 0.1038738 
kable(coef(h1_sex_mod_dataset_varying)$dataset[, ,"f0"], caption = "Estimated association between f0 and outcome by dataset")
Estimated association between f0 and outcome by dataset
Estimate Est.Error Q2.5 Q97.5
2 0.1664338 0.1095638 -0.0368466 0.3912092
3 0.0060973 0.1109487 -0.2197120 0.2172270
4 0.0183517 0.1379504 -0.2793900 0.2752177
7 0.0456925 0.1420845 -0.2490861 0.3217105
8 -0.0038451 0.1582957 -0.3755189 0.2684981
9 0.0060925 0.1247970 -0.2605025 0.2344679
10 0.1615680 0.1726996 -0.1191796 0.5705388

Hypothesis 4: SOI behavior

Participants with lower voice pitch will report having a more unrestricted sociosexual behavior.

Interpretation: The data are are consistent with substantial linear negative relationship between f0 and unrestricted sociosexuality. The 89% HDI for f0 falls entirely outside the ROPE, but the HDI for pf falls almost entirely within it. This evidence is consistent with a non-negligible association, where people with deeper voices have a less restricted sociosexuality (after adjusting for age and gender). Further research is needed to verify if the association with pf is truly negligible in size.

Visually diagnose non-linearity

ggplot(vcs, aes(f0, behavior)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(pf, behavior)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex") 

ggplot(vcs, aes(age*10, behavior)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Decisions upon visual diagnosis: We will fit splines for f0, pf, and age. It seems likely that the nonlinearities in f0 and pf will not both be apparent after adjusting for one another. There is no evidence in the plots that there is an interaction with sex, except for age.

Models

h1_simple <- brm(behavior ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/behavior/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 


h1_after_diagnosis <- brm(behavior ~  s(f0) + s(pf) + sex + s(age, by = sex) + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, cores = chains, chains = chains,
          prior = priors, 
          control = control, save_mevars = TRUE,
          file = "models/behavior/h1_nonlinear") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
conditional_smooths(h1_after_diagnosis)

compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
##                    elpd_diff se_diff
## h1_after_diagnosis  0.0       0.0   
## h1_simple          -9.2       5.0
h1_sex_mod_dataset_varying <- brm(behavior ~ f0 + sex + s(age, by = sex) + 
                                    (1 + f0 + sex || dataset), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/behavior/h1_sex_mod_dataset_varying2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_after_diagnosis_melsm <- brm(
  bf(behavior ~  f0 + sex + s(age, by = sex) + (1 | dataset),
     sigma ~ f0 + sex + s(age, by = sex) + (1 | dataset)), data = vcs,
          prior = priors,
          iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
          control = control, save_mevars = TRUE,
          file = "models/behavior/h1_after_diagnosis_melsm") %>%
  add_criterion("loo") %>%
  add_criterion("bayes_R2") %>%
  add_criterion("loo_R2")

compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying, h1_after_diagnosis_melsm)
compare_models
##                            elpd_diff se_diff
## h1_after_diagnosis_melsm     0.0       0.0  
## h1_sex_mod_dataset_varying -74.7      15.0  
## h1_after_diagnosis         -75.1      14.9  
## h1_simple                  -84.3      15.8
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying, h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>% 
  map(~ { x <- as_tibble(bayes_R2(.))
          x$loo_R2 <- loo_R2(.)
          x
          }) %>%
  bind_rows(.id="model") %>% 
  arrange(desc(loo_R2))

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)
  behavior
Predictors Estimates CI (95%)
Intercept -0.84 -1.18;-0.49
Voice pitch (f0) -0.20 -0.32;-0.09
Formant position (Pf) -0.03 -0.15;0.10
Sex [male] -0.24 -0.39;-0.09
Age 0.33 0.24;0.42
Random Effects
σ2 0.10
Ï„00 0.90
ICC 0.10
N dataset 9
Observations 1996
Marginal R2 / Conditional R2 0.154 / 0.154 
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter H0 inside ROPE 89% HDI
Intercept Rejected 0% [-1.10;-0.57]
Voice pitch (f0) Rejected 0% [-0.30;-0.11]
Formant position (Pf) Undecided 91% [-0.13; 0.07]
Sex [male] Rejected 0% [-0.36;-0.12]
Age Rejected 0% [ 0.26; 0.40] 
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Model details

summary(h1_simple)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: behavior ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1996) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.36      0.12     0.21     0.67 1.00     2139     3425
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.84      0.17    -1.18    -0.49 1.00     2843     4344
## f0           -0.20      0.06    -0.32    -0.09 1.00     5532     5261
## pf           -0.03      0.06    -0.15     0.10 1.00     6272     5436
## sex1         -0.24      0.08    -0.39    -0.09 1.00     4784     4804
## age           0.33      0.05     0.24     0.42 1.00     6533     5055
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.92      0.01     0.89     0.95 1.00     7919     5519
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis)
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: behavior ~ s(f0) + s(pf) + sex + s(age, by = sex) + (1 | dataset) 
##    Data: vcs (Number of observations: 1996) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Smooth Terms: 
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sf0_1)               0.76      0.74     0.03     2.71 1.00     1880     3952
## sds(spf_1)               0.58      0.55     0.02     1.98 1.00     3006     4303
## sds(sagesexfemale_1)     1.83      0.90     0.54     3.95 1.00     3379     3912
## sds(sagesexmale_1)       1.20      0.78     0.21     3.20 1.00     3019     3731
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.35      0.12     0.20     0.64 1.00     2630     4222
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           -0.03      0.12    -0.28     0.22 1.00     1671     2747
## sex1                -0.24      0.10    -0.44    -0.04 1.00     6988     5189
## sf0_1               -1.44      1.40    -4.00     2.04 1.00     4390     3647
## spf_1               -0.05      1.25    -2.42     2.85 1.00     5874     4280
## sage:sexfemale_1     1.70      2.27    -3.16     5.89 1.00     4914     5346
## sage:sexmale_1       0.62      1.69    -3.23     3.65 1.00     4798     4744
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.92      0.01     0.89     0.94 1.00    11917     6138
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(h1_simple)

conditional_smooths(h1_after_diagnosis)

summary(h1_sex_mod_dataset_varying)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: behavior ~ f0 + sex + s(age, by = sex) + (1 + f0 + sex || dataset) 
##    Data: vcs (Number of observations: 1996) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Smooth Terms: 
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sagesexfemale_1)     1.79      0.90     0.51     3.95 1.00     4849     5303
## sds(sagesexmale_1)       1.21      0.79     0.20     3.17 1.00     4511     5097
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.35      0.13     0.18     0.66 1.00     5788     7785
## sd(f0)            0.09      0.08     0.00     0.28 1.00     4223     7043
## sd(sex1)          0.10      0.10     0.00     0.38 1.00     3743     5264
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           -0.04      0.13    -0.31     0.23 1.00     5186     6791
## f0                  -0.20      0.07    -0.34    -0.06 1.00    13673     8382
## sex1                -0.22      0.09    -0.41    -0.03 1.00     9046     6742
## sage:sexfemale_1     1.73      2.24    -3.06     5.85 1.00    10938     8981
## sage:sexmale_1       0.67      1.67    -2.98     3.78 1.00     9976     8620
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.92      0.01     0.89     0.95 1.00    24933     8698
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis_melsm)
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: behavior ~ f0 + sex + s(age, by = sex) + (1 | dataset) 
##          sigma ~ f0 + sex + s(age, by = sex) + (1 | dataset)
##    Data: vcs (Number of observations: 1996) 
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup samples = 12000
## 
## Smooth Terms: 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sagesexfemale_1)           1.46      0.82     0.34     3.47 1.00     5871     6439
## sds(sagesexmale_1)             1.04      0.70     0.18     2.85 1.00     5497     5796
## sds(sigma_sagesexfemale_1)     0.56      0.51     0.02     1.88 1.00     4447     5846
## sds(sigma_sagesexmale_1)       0.76      0.53     0.10     2.10 1.00     5619     5142
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.36      0.12     0.21     0.67 1.00     4431     6769
## sd(sigma_Intercept)     0.27      0.09     0.15     0.50 1.00     4616     6540
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept                 -0.02      0.13    -0.29     0.23 1.00     3639     5512
## sigma_Intercept           -0.15      0.10    -0.35     0.05 1.00     3955     5418
## f0                        -0.18      0.05    -0.29    -0.08 1.00    11943     8852
## sex1                      -0.19      0.06    -0.30    -0.08 1.00    12981     8786
## sigma_f0                  -0.01      0.05    -0.09     0.08 1.00    13091     8714
## sigma_sex1                 0.05      0.05    -0.05     0.15 1.00    12586     9044
## sage:sexfemale_1           1.78      2.09    -2.73     5.64 1.00    10448     8689
## sage:sexmale_1             0.53      1.60    -2.97     3.37 1.00     9531     8351
## sigma_sage:sexfemale_1     0.18      1.46    -3.11     3.13 1.00     6791     5141
## sigma_sage:sexmale_1       0.47      1.46    -2.32     3.77 1.00     8548     7422
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
Explained variance, regular and leave-one-out.
model Estimate Est.Error Q2.5 Q97.5 loo_R2
h1_after_diagnosis_melsm 0.1625449 0.0131309 0.1374574 0.1888755 0.1539330
h1_sex_mod_dataset_varying 0.1676392 0.0138851 0.1403527 0.1949088 0.1532208
h1_after_diagnosis 0.1682773 0.0142038 0.1408193 0.1965213 0.1527536
h1_simple 0.1542048 0.0138106 0.1269713 0.1812360 0.1451403 
kable(coef(h1_sex_mod_dataset_varying)$dataset[, ,"f0"], caption = "Estimated association between f0 and outcome by dataset")
Estimated association between f0 and outcome by dataset
Estimate Est.Error Q2.5 Q97.5
1 -0.2042854 0.0788527 -0.3624871 -0.0472227
2 -0.2490503 0.0929913 -0.4585575 -0.0898764
3 -0.1656330 0.0927578 -0.3309920 0.0412035
4 -0.1530315 0.1044902 -0.3264683 0.0901174
5 -0.2028399 0.0898606 -0.3849846 -0.0220809
6 -0.1937525 0.0996324 -0.3902962 0.0127610
7 -0.2306656 0.1075203 -0.4808581 -0.0413506
8 -0.1662645 0.1085892 -0.3568731 0.0883545
11 -0.2264906 0.1189240 -0.5077337 -0.0120300 
equivalence_test(h1_after_diagnosis_melsm, range = c(-0.1, 0.1))
## # Test for Practical Equivalence
## 
##   ROPE: [-0.10 0.10]
## 
## Parameter              |        H0 | inside ROPE |       89% HDI
## ----------------------------------------------------------------
## Intercept              | Undecided |     67.86 % | [-0.22  0.18]
## sigma_Intercept        | Undecided |     26.54 % | [-0.31 -0.00]
## f0                     | Undecided |      0.04 % | [-0.27 -0.10]
## sex1                   |  Rejected |      0.00 % | [-0.29 -0.10]
## sigma_f0               |  Accepted |    100.00 % | [-0.08  0.07]
## sigma_sex1             | Undecided |     89.17 % | [-0.03  0.13]
## sage.sexfemale_1       | Undecided |      2.35 % | [-1.52  5.01]
## sage.sexmale_1         | Undecided |      4.75 % | [-2.08  2.92]
## sigma_sage.sexfemale_1 | Undecided |      8.28 % | [-2.19  2.16]
## sigma_sage.sexmale_1   | Undecided |      7.80 % | [-1.82  2.73]

Summary of results

Models

dominance <- readRDS("models/dominance/h1_simple.rds")
agree <- readRDS("models/agree/h1_simple2.rds")
extra <- readRDS("models/extra/h1_simple.rds")
soi <- readRDS("models/behavior/h1_simple.rds")
theme_set(theme_cowplot())

modify <- function(outcome) {
  list(scale_y_discrete(outcome, breaks = names(effect_labels), labels = effect_labels),
  guides(fill = "none"),
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0, 0.5, 0, 0.5), "cm"))
)
}
remove_x_axis <- list(
  scale_x_continuous("", limits = c(-0.5, 0.4)),
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        axis.line.x = element_blank())
)

p_dominance <- plot(equivalence_test(dominance, parameters = c("b_f0", "b_pf", "b_sex1"))) + modify("Dominance") + remove_x_axis
p_agree <- plot(equivalence_test(agree, parameters = c("b_f0", "b_pf", "b_sex1"))) +
  modify("Agreeableness") + remove_x_axis
p_extra <- plot(equivalence_test(extra, parameters = c("b_f0", "b_pf", "b_sex1")))  +
  modify("Extraversion") + remove_x_axis
p_soi <- plot(equivalence_test(soi, parameters = c("b_f0", "b_pf", "b_sex1"))) +
  modify("SOI behavior") + xlim(c(-0.5, 0.4))
legend <- get_legend(
  p_dominance + 
    guides(fill = guide_legend(nrow = 1)) +
    theme(legend.position = "bottom", legend.justification = "center")
)
## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).
neg_spacer <- -0.17
aligned_plots <- align_plots(p_dominance, p_extra, p_agree, p_soi, align="hv")
## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).
plot_grid(ncol = 1,  rel_heights = c(1, neg_spacer, 0.8, neg_spacer, 1, neg_spacer, 1.1, .15),
          aligned_plots[[1]], NULL, aligned_plots[[2]], NULL, aligned_plots[[3]],NULL, aligned_plots[[4]], legend)

ggsave("preregistered_rope.pdf", width = 8, height = 10)

Tabular

get_coefs <- function(model) {
model_summary <- summary(model)
nonvarying <- model_summary$fixed %>% 
  as.data.frame() %>% 
  rownames_to_column("term") %>% 
  mutate(group = "non-varying",
         outcome = as.character(model_summary$formula[[1]][[2]]))
varying <- model_summary$random %>% 
  map(~ rownames_to_column(as.data.frame(.), "term")) %>% 
  bind_rows(.id = "group") %>% 
  left_join(model_summary$ngrps %>% 
    as_tibble() %>% 
    gather(group, n), by = "group") %>% 
  mutate(#group = paste0(group, " (n=", n, ")"),
         outcome = as.character(model_summary$formula[[1]][[2]]))

two_digits <- function(x) { sprintf("%.2f", x) }
coefs <- bind_rows(nonvarying, 
                   varying)
coefs <- coefs %>% select(group, outcome, term, Estimate, `l-95% CI`, `u-95% CI`)
coefs$Estimate = two_digits(coefs$Estimate)
coefs$`l-95% CI` = two_digits(coefs$`l-95% CI`)
coefs$`u-95% CI` = two_digits(coefs$`u-95% CI`)
coefs
}

coefs <- bind_rows(get_coefs(dominance), 
                   get_coefs(extra), get_coefs(agree), get_coefs(soi))
## Warning: There were 3 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
wide_coefs <- coefs %>% 
  mutate(Estimate = paste0(Estimate, " [", `l-95% CI`, ";", `u-95% CI`, "]")) %>% 
    select(group, outcome, term, Estimate) %>% 
  mutate(term = recode(term, "meageage_se" = "Age(±SE)", "age" = "Age(±SE)"),
         outcome = fct_inorder(outcome),
         term = fct_inorder(term)) %>% 
  spread(outcome, Estimate) %>% 
  arrange(desc(group), term) %>% 
  mutate(group = fct_inorder(group))

outcome_names <- c("Extraversion" = "extra", "Dominance" = "dominance",
                   "Agreeableness" = "agree", "SOI behavior" = "behavior")

library(kableExtra)
wide_coefs %>% 
  mutate(term = str_replace_all(term, variable_labels)) %>% 
  rename(!!!syms(outcome_names), Term = term) %>% 
  select(-group) %>% 
  kable(caption = "Exploratory analysis results") %>%
  pack_rows(index = table(wide_coefs$group)) %>% 
  column_spec(column = 1, width = "3.3cm") %>% 
  column_spec(column = 2:(ncol(wide_coefs)-1), width = "3cm") %>% 
  add_header_above(c(" " = 1, "Estimated effect on each outcome [95% CI]" = ncol(wide_coefs)-2)) %>% 
  footnote("Estimated associations between outcomes, for which we had made preregistered predictions.")
Exploratory analysis results
Estimated effect on each outcome [95% CI]
Term Dominance Extraversion Agreeableness SOI behavior
non-varying
Intercept -0.07 [-0.75;0.60] 0.17 [-0.22;0.57] 0.08 [-0.38;0.56] -0.84 [-1.18;-0.49]
Voice pitch (f0) -0.27 [-0.45;-0.09] -0.23 [-0.38;-0.08] 0.02 [-0.13;0.17] -0.20 [-0.32;-0.09]
Formant position (Pf) -0.02 [-0.20;0.16] 0.08 [-0.05;0.21] 0.04 [-0.09;0.17] -0.03 [-0.15;0.10]
Sex [male] -0.31 [-0.53;-0.09] -0.30 [-0.48;-0.11] -0.01 [-0.18;0.17] -0.24 [-0.39;-0.09]
Age(±SE) 0.04 [-0.08;0.16] -0.10 [-0.22;0.01] -0.06 [-0.18;0.05] 0.33 [0.24;0.42]
dataset
sd(Intercept) 0.51 [0.16;1.52] 0.35 [0.16;0.74] 0.47 [0.24;0.98] 0.36 [0.21;0.67]
Note:
Estimated associations between outcomes, for which we had made preregistered predictions.

Scatterplots

dominance <- ggplot(vcs, aes(f0, dominance)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("") +
  ylab("Dominance") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))

extra <- ggplot(vcs, aes(f0, extra)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("") +
  ylab("Extraversion") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))

agree <- ggplot(vcs, aes(f0, agree)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("Voice pitch (f0)") +
  ylab("Agreeableness") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))


behavior <- ggplot(vcs, aes(f0, behavior)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("Voice pitch (f0)") +
  ylab("SOI Behavior") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))


plot_grid(ncol = 2,  dominance, extra, agree, behavior, rel_heights = c(1, 1, 1, 1),
          align = "hv")
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing missing values (geom_point).
## Warning: Removed 797 rows containing non-finite values (stat_smooth).
## Warning: Removed 797 rows containing missing values (geom_point).
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).
## Warning: Removed 234 rows containing non-finite values (stat_smooth).
## Warning: Removed 234 rows containing missing values (geom_point).

ggsave("scatterplots_prereg.pdf", width = 8, height = 8)